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ABSTRACT 

Aims. We present a spectroscopic study of a sample of 287 planetary nebulas (PNs) around the brightest cluster galaxy (BCG) M87 in 
Virgo A, of which 211 are located between 40 kpc and 150 kpc from the galaxy centre. With these data we can distinguish the stellar 
halo from the co-spatial intraciuster light (ICL) and study both components separately. 

Methods. We obtained PN velocities with a high resolution ELAMESA^LT survey targeting eight fields in a total area of ~ 0.4 deg“. 
We identified PNs from their narrow and symmetric redshifted /15007A [GUI] emission line, the presence of the second /14959A [GUI] 
emission line, and the absence of significant continuum. We implement a robust technique to measure the halo velocity dispersion 
from the projected phase-space to identify PNs associated with the M87 halo and ICL. Using photometric magnitudes, we construct 
PN luminosity functions (PNLEs), which are complete down to m5oo7=28.8. 

Results. The velocity distribution of the spectroscopically confirmed PNs is bimodal, containing a narrow component centred on the 
systemic velocity of the BCG and an off-centred broader component, which we identify as halo and ICL, respectively. We find that 
243 PNs are part of the velocity distribution of the M87 halo, while the remaining subsample of 44 PNs are intraciuster PNs (ICPNs). 
Halo and ICPNs have different spatial distributions: the number density of halo PNs follow the galaxy’s surface brightness profile, 
whereas the ICPNs are characterised by a shallower power-law profile, Lcl with y in the range [-0.34, -0.04]. No evidence is 
found for an asymmetry in the halo and ICPN density distributions when the NW and SE fields are studied separately. A study of the 
composite PN number density profile confirms the superposition of different PN populations associated with the M87 halo and the 
ICL, characterised by different PN specific numbers a. We derive tthaio = L06 x 10“** Np^L^^'j^^j and q-icl = 2.72 x I0“* * Np^L^'i^^i, 
respectively. The M87 halo PNLE has fewer bright PNs and a steeper slope towards faint magnitudes than the ICPNLE, and both are 
steeper than the standard PNLE for the M31 bulge. Moreover, the ICPNLE has a dip at ~ 1-1.5 mag fainter than the bright cut-off, 
reminiscent of the PNLEs of systems with extended star formation history, such as M33 or the Magellanic clouds. 

Conclusions. The BCG halo of M87 and the Virgo ICL are dynamically distinct components with different density profiles and 
velocity distributions. Moreover, the different a-parameter values and PNLE shapes of the halo and ICL indicate distinct parent stellar 
populations, consistent with the existence of a gradient towards bluer colours at large radii. These results reflect the hierarchical 
build-up of the Virgo cluster. 
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1. Introduction 

Galaxy halos are faint stellar components made of stars grav¬ 
itationally bound to the individual galaxies. In galaxy clusters 
these halos may be surrounded by intraciuster stars. The exis¬ 
tence of a diffuse population o f intergalactic stars was first pro¬ 
posed bv lZwick'^(ll937lll952ll . As a consequence of its low sur¬ 
face brightness, it was only with the advent of CCD photometry 
that this diffuse stellar component could be studied in a quantita¬ 
tive way, thus becoming a topic of interest for observational and 
theoretical studies. 

The formation of intraciuster light (ICL) and of the extended 
halos around the brightest cluster galaxies (BCG) is closely re¬ 
lated to the morphological transformation of galaxies in clus¬ 
ters. Two of the main physical processes describing the gravita- 
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tional interaction between galaxies during cluster formation and 


Merritt 19851; 

Taffoni et alJ l2003l iBovlan-Kolchin et al. 

2008; 

De Lucia et al. 

2 OI 2 I) and tidal stripping ( 

Gallagher & Ostrikei 

1972 

: iMoore et alJ Il996l iGress 8c WestI 

19981: Willmanetal. 


tance, and orbit, these processes determine the fate of a clus¬ 
ter galaxy. Dynamical friction is the primary mechanism drag¬ 
ging a massive satellite towards the host halos centre, where 
it merges with the BCG. On the other hand, tidal forces strip 
stars from satellite galaxies which end up orbiting the clus- 
ter as unbou nd ICL dGnedinl 200^ iMurante et alJl20()4l I2007I : 
lMihosll2004t iRudick et al.ll2006 ). In the ICL and the outer re¬ 
gions of BCGs where the dynamical timescales are long, fos¬ 
sil recor ds of accretion events can be preserved over extended 
perio ds dWillman et alJlTOOil iRudick et al.l |200^ ICooper et al.l 
l2014l) . Hence, the study of the luminosity, distribution, and kine¬ 
matics of galaxy halos and ICL may provide information on the 
evolution of galaxies and their host clusters. 
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In the last ten years, 
clusters has shed light 


the analysis of simulated galaxy 
upon the nature and origin of 


BCG dNanolitano et alJ l2003l iMurante et al. 

2004 2007t 

Willman et al.ll2004]; iRudick et al.l 20061 20091; 

Puchwein et all 

2 OIOI; Lanorte et alJ l2013l IContini et al.l 20141) 

In the frame- 


work of co smologica l hydro dynamic al simulations. iDolag et al.l 
( I 2 OIOI) and ICui et al.1 (1201 4h separated stars bound to the clus¬ 
ter potential from those bound to the BCG by adopting a dy¬ 
namical dehnition of the main galaxy halo and the diffuse light. 
Tagging particles as galaxy or intracluster component based on 
their different yelocity distributions, they identihed two distinct 
stellar populations in terms of kinematics, spatial distribution, 
and physical properties like age and metallicity. Other studies 
do not adopt this kind of dynamical dehnition of the two com¬ 
ponents and treat the BCG and ICL as a single system at the 
centre of the cluster, consisting of all star s that are not bound 
to any satellite sub halos in the cluster (e.g. lMurante et al1l2004t 
ICooner et al .11201^ . 

Galaxy halos, diffuse ICL, and their connection with galaxy 


ies. Deep imag 

rng of indiyidual objects (iBernstein et al. 

1995 

Gonzalez et al.l 

120001 Feldmeier et al. 

2004 iMihos et al. 

2005 

Krick & Bernsteinl2007 ;lRudick et al. 

20101) shows that the faint 


ICL around BCGs often has irregular morphology, consistent 
with predictions from simul ations. From stack i ng images for 
a larg e number of objects (IZibetti et al.l l200.5^ ID’ Souza et al.l 
|20 i 4 . ayerage photometric properties were obtained, showing 
that the ICL extends to many 100s of kpcs from the cluster cen¬ 
tre. In the Virgo cluster. iKormendy et al.l(l2009h analysed a sam¬ 
ple of ellipticals and spheroidal galaxies and studied their halos 
through the light profiles. Comparing their results with the earli¬ 
est simulations of lyan Albadal (Il982h they argued that light pro- 
hles with a large Sersic index (n > 4) are common in many giant 
ellipticals whose origin can be associated with merger processes. 

Kinematics and stellar population parameters haye only been 
measured in a small number of BCG halo s, such as (i) NGC 4889 
in the Coma cluster (ICoccato et al.ll2010l) . which shows a change 
of stellar popula tion at large radii, (ii) the central galaxy NGC 
3311 in Hydra I (IVentimiglia et al.l 2010t[Arnaboldi et al.ll2012]) , 
where the kinematics as well as the morphology of the outer 
halo signal on-going accretion eyents, and (hi) N GC 6166 in 
Abell 2199 (iKelson et al.ll2002t iBender et al.ll2014l) , whose ye¬ 
locity dispersion prohle blends smoothly into the cluster. In the 
Coma cluster, the ICL kinematics sugg est an on-going merger of 
two cluster cores (iGerhard et al.ll20()% . All these studies aimed 
to understand the role of tidal disruption and merger eyents as 
the main processes inyolyed in the formation and eyolution of 
central cluster galaxies and the ICL. 

In the nearest clusters, single stars can be used to study 
the stellar populations asso ciated with the outer halos and the 
diffuse stellar component (|Fergu^ro et^.| 1998t iDurrell et al.l 
l2002t [williams et al.ll2007[ lYan et 12008 ). Globular clusters 
(GCs) haye been used to obtain kinematic i nformation in the 
outer regions of near b y early-type galaxi es (Cote_e^^ 20011 
Schuberth et al. 12010: IStrader et al.l 1201 ll: iRomanowsky et akl 


2012 ! iPota et al] l2013h . Planetary Nebulas (PNs) haye been 
targeted in seyeral sury eys aimed to tra c e the light and mo 


tions in 


Peng et al.l 


galactic 1 

u [2001 


'eys _ _ _ 

a lactic h alos (IHui et al.[ Il993t jMendez_et^ 

' 200T 


2001 


201c; 


_ _ jCoccato et al.l I2009t iMcNeil et ak. _ 

McNeil-Moylan et al. 2012t Cortesi et aiT2013h. the Virgo 
cluster IC component dA rnaboldi et al.l 1 1 9961 200^ 20o4 


Agueiri et al.ll2005l: Doherty et al.ll2009HCastro-Rodrfguez et alJ 

2009Hf^ngobardi et al.ll2013h . and the Hydra I and Coma clus¬ 


ters, out to 50-100 Mpc distance (IVentimiglia et al.ll2010l 120111 
iGerhard et al.l l200.5ll . It was found that the obseryed proper¬ 
ties of the PN population, such as the a parameter that quan- 
tihes the stellar luminosity associated with a detected PN, and 
the PN luminosity function (PNLF) correlate with the age, 
colour, and metallicity of the parent stellar population (Hui et all 


colour, and metallicity of the parent stellar population dHui et al. 

1993HCiardullo et al.l2004ICiardullol2010tlBuzzoni et al.l2006 


Longobardi et al.N2013l) . Thus PNs can be used to trace these 

physical quantities of their parent stellar populations at surface 
brightnesses too faint for other techniques. 

The giant elliptical galaxy M87 has one of the oldes t stel¬ 
lar populations in the local Uniyerse (iLiu et al.l l2005h . and 
a stellar halo containin g 70% of the galaxy light down to 
fly -27.0 mag arcsec^^ dlformendy et al.ll2009ll. It is close to 
the ce ntre of sub-cluster A in the Virgo cluster (iBinggeli et al.l 
1 19871) . the nearest galaxy cluster, and it is expected to haye 
transformed oyer larger tim escales because of galaxy mergers 
(iDe Lucia & Blaiz^l2007l) . Deep imaging (iMihos etalJ 120^051 
iJanowiecki et al.ll201ol) has reyealed a complex network of faint, 
extended tidal features around M87, suggesting that it is not 
completely in equilibrium. Thus, M87 and the surrounding Virgo 
cluster core are prime targets to address the formation and 
eyolution of galaxy clusters, ICL, and BCGs. Indeed, M87 is 
the su bject of many dynamical studies with X-ray me asure- 
ments dNulsen & Bobrin^ 1995; Churazoy et alJ TOlOh . inte- 


matics 


ics (|Cote et al.||200 1l:IStrader et alj2()lTF RomanowskY_et^ 


ments (INulsen & Dohrmgerlliqqpt IChurazoy et alJI/t)t(JI) . inte¬ 
grated stellar kinematics lMumhy-5!a L 201 iL 20141). GC kine- 


20.12 ; |Zhu et al.l 20141) a nd PN kinematics ( Arnaboldi et al 
2004 : liberty et al. l2009l) . to estimate its mass and deriye the 
dark m atter distribution. Using PN kinematics, iDoherty et al.l 
(I 2 OO 9 I) identihed M87 halo and IC PNs and showed the coex¬ 
istence, at radii > 60 kpc, of a stellar halo bound to the galaxy 
potential and a surrounding unbound Virgo ICL. 

In this work, we report the results of a wide and high res¬ 
olution spectroscopic suryey coyering the outer regions of M87 
out to a distance of 150 kpc from the galaxy centre. The aim 
of this project is to inyestigate the halo-ICL dichotomy, making 
use of a large spectroscopic sample of P Ns (approximately 1 5 
times larger than the preyious sample of iDoherty et al.ll2009l) . 
The paper is structured as follows; in Sect. |2] we describe the 
spectroscopic suryey together with the data reduction procedures 
and the classihcation of PN spectra. In Sect. [3 we study the 
PN phase-space distribution and dynamically separate the halo 
and IC populations. Spatial density distributions are deriyed in 
Sect.|4] and in Sect.|5]we present the properties of the halo and 
IC PN populations in terms of their a-parameters and the mor¬ 
phology of their PNLFs. Finally, we discuss our results in Sect.|6] 
and giye our conclusions in Sect. [T] 

We adopt a _ dist ance modulus of 30 .8 for M87 

(ICiardullo et al.l l2002t iLongobardi et al.l 12013l) . implying a 
physical scale of 73 pc arcsec^*. 


2. The FLAMES M87 PN survey 

2.1. Photometric sample 

The photometric candidates targeted by our s pectroscopic sur¬ 
yey c ome from an earlier imaging suryey (ILongobardi et al.l 
I 2 OI 3 I) . coyering a 0.43 deg^ region centred on M87. Images 
were taken through a narrowband hlter centred on the redshifted 
[OIII]/15007A emission line at the Virgo cluster distance (on- 
band image), and through a broadband V-hlter (off-hand image). 
Because of their bright [OIII]/15007A emission, extragalactic 
PNs can be identified as unresolyed emission sources with pos- 
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itive flux on the on-band [OIII] image and no detection on the 
off-hand image. 

We obtained spectra in two observing campaigns. For the 
first spectroscopic campaign, we selected emission line candi¬ 
dates as objects with positive flux on th e colour tOnil-V band 
imag43 (hereafter difference method, see lFeldmeier et all (l2003h 
for more details). The visual catalogue extracted on the basis of 
the difference method consisted of 1074 objects, which covered 
a magnitude range 23 < ^ 29.8, and is statistically com¬ 

plete down to OT 5007 - 28.® 

For the second spectroscopic campaign, we carried 
out a more st r ingent selection procedure described in 
iLongobardi et alJ (l2Q13h . In this procedure, the PN candi¬ 
dates were selected using automatic selection criteria, based 
on the distribution of the detected sources in the colour- 
magnitude (CM) diagram and the pro perties of their point- 
spread function (PSF; fo r more details see lArnaboldi et al.l2002t 
ILongobardi et alj|2013h . This automatic catalogue is complete 
within the magnitude range 26.3 < m 5007 < 28.4. 

The combined total input sample for the spectroscopic sur¬ 
vey (visual catalogue plus automatic catalogue) consisted of 
1484 emission line candidates. 

2.2. Observations and data reduction 

The spectra were acquired in service mode with the 
FLAMES spectrograph on the VLT-UT2 telescope, in the 
GIRAFFEh-MEDUSA mode. This observing mode allows for 
up to 132 separate fibres that can be allocated to targets in one 
plate configuration, covering a circular area of 20' diameter. The 
total emission line sample was observed in two observing runs 
(24h, 088.B-0288(A); llh, 093.B-0066(A); PI; M.Arnaboldi), 
which were characterised by clear conditions and seeing better 
than 0.9". We used the high resolution grism HR08, covering a 
wavelength range of ~ 250 A centred on 5048 A with a spec¬ 
tral resolution of R = 22500. With this setup, the instrumental 
broadening of the arc lines has a EWHM of 17 km s * and the 
stati stical error on the wavelength measurements is 150 m s“' 
(see iRover et al.ll2002h . We refer to Sect. 12.51 for discussion of 
the velocity accuracy estimated from repeat observations of the 
same emission line candidates. 

Because the [OIII] emission lines from PNs are only a few 
km s * wide, high resolution spectra are also desirable to reduce 
the sky contamination, making the ELAMES spectrograph the 
ideal instrument for LOS velocity measurements of extragalactic 
PNs. 

Eor our first spectroscopic campaign, the visual catalogue 
was divided into a bright (msoor < 27.2) and a normal (msoor > 
27.2) sample. The ELAMES plate configurations and exposure 
times were then optimised to reach the maximum number of 
fibres allocated, as well as optimal signal-to-noise ratio (S/N) 
for both samples. Based on the ELAMES/GIRAEEE/HR08 
configuration, for a lO'^ sec exposure, the monochromatic 
[OIII]/15007A emission of a PN with [OIII] flux of P 5007 = 
1.0xl0“*^ergcm“^s“* (m5oo7 = 28.8) is detected with S/N= 10 
per resolution element (0.29 A). 


* Before the subtraction the continuum off-band image was scaled to 
the on-band image by a multiplicative scaling factor found measuring 
fluxes from several bright, isolated stars on both images. 

^ The transformation between the AB and 5007 magnitudes for 
the photometric narrow band filter is given by 025007 = + 2.49 

dLongobardi et alJ2013h . 



Right ascension 



188.2 188.0 187.8 187.6 187.4 187.2 

R.A. (degree) 


|1g.l. Top panel: Image of the core of the Virgo cluster dMihos et all 
120051) with the positions of the FLAMES fields studied in this 
work (cyan circ l es) an d in previous surveys (dot ted black circles) by 
lArnaboldi et all d2004t) and lPohertv et all d2009l) . Black squares rep¬ 
resent the Suprime-Cam sur vey fields used for the extraction of the 
photometric PN candidates dLongobardi et all l2013h . Bottom panel: 
Schematic zoom-in of the surveyed area (black rectangles). FLAMES 
pointings (cyan circles) and PN candidat es from the photometric survey 
(black asterisks: ILongobardi et all |2013[) ar e plotted over th e isophotes 
of the M87 halo surface brightness dKormendv et al1l2009t) . The green 
cross depicts M87’s centre. North is up. East is to the left. 


In total, we defined 12 FLAMES plate configurations, la¬ 
belled as M87SUB1 Bright F01-F03, M87SUB1 Norm F01-F03 
and M87SUB2 Bright F01-F03, M87SUB2 Norm F01-F03 for 
the NW (SUBl) and SE (SUB2) Suprime-Cam fields, respec¬ 
tively. The layout of the FLAMES pointings on the sky, together 
with the coverage of the photometric Suprime-Cam survey, is 
shown in Fig.[T] 

In the second spectroscopic campaign, we added two 
FLAMES plate configurations covering the very central region 
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Table 1. Flames configuration and exposure times. Column 1: ID of FLAMES plate. Columns 2 & 3: FLAMES plate; right ascension and 
declination. Column 4; Total exposure time. Column 5: Single exposure time x number of exposures. 


FLAMES Conf 

RA 

h:m:s 

DEC 

degrees 

Total Exposure time 

(s) 

Single Exposure time 

X # of exposures 
(s) 

M87SUB1 Bright FOl 

12:30:58.015 

-h12:39:39.28 

2700 

1350 

X 2 

M87SUB1 Bright F02 

12:30:05.332 

-h12:34:00.08 

2700 

1350 

X 2 

M87SUB1 Bright F03 

12:29:59.574 

+ 12:30:34.49 

2700 

1350 

X 2 

M87SUB1 Norm FOl 

12:30:58.015 

+ 12:39:39.28 

8100 

2700 

X 3 

M87SUB1 NormF02 

12:30:05.332 

+ 12:34:00.08 

8400 

2800 

x3 

M87SUB1 NormF03 

12:29:59.574 

+ 12:30:34.49 

8400 

2800 

X 3 

M87SUB2 Bright FOl 

12:31:25.426 

+ 11:57:31.43 

2700 

1350 

X 2 

M87SUB2 Bright F02 

12:31:44.477 

+ 12:12:01.55 

2700 

1350 

X 2 

M87SUB2 Bright F03 

12:30:45.170 

+ 12:02:27.20 

2700 

1350 

X 2 

M87SUB2 Norm FOl 

12:31:25.426 

+ 11:57:31.43 

8400 

2800 

X 3 

M87SUB2 Norm F02 

12:31:44.477 

+ 12:12:01.55 

8400 

2800 

X 3 

M87SUB2 Norm F03 

12:30:45.170 

+ 12:02:27.20 

8400 

2800 

X 3 

M87SUB FEDGE 

12:29:57.854 

+ 12:41:20.90 

6800 

1700 

x4 

M87SUB FCEN 

12:30:50.019 

+ 12:20:48.88 

6800 

1700 

x4 


of M87 and the NW edge of the Suprime-Cam imaging survey, 
in addition to completing the observations of the 12 FLAMES 
configurations from the first campaign. These two additional 
FLAMES configurations are shown in Fig. [1] with the labels 
FCEN and FEDGE, respectively. Table [1] provides an overview 
of the FLAMES field configurations and the total exposure 
times. 

We reduced the spectroscopic data using the GIRAFFE 
pipelin^B- The reduction procedure included bias subtraction, 
flat-fielding, identification of the fibre locations on the CCD, 
geometric distortion correction, wavelength calibration and ex¬ 
traction of the one-dimensional spectra. The calibrated one¬ 
dimensional spectra were then corrected to the heliocentric ve¬ 
locity using the IRAFQ task dopcor. Finally, we combined spec¬ 
tra from single exposures using the IRAF task scombine to get 
the targeted S/N for each spectrum (see above). 

2.3. Spectroscopic success rates 

We define the nominal success rate as the ratio between the 
number of spectra with a detected emission line and the num¬ 
ber of fibres allocated for a given FLAMES plate configuration. 
For our observations, it varies from field to field, in a range of 
values between ~ 20% and ~ 60%. These values are similar 
to those obtained in the spectro scopic follow - up of GCs in the 
outer halo of M87 (see Sect. 3 in IStrader et al.l(l201 ll) "Pl The low 
success rates for some of the fields (mostly M87SUB2 fields) 
are caused by guide star proper motions, which were not cor¬ 
rectly accounted for in the FLAMES astrometry input file for 
the fibre allocation. We were also able to estimate the fraction of 
fibre-object misalignments, from the repeat observations of the 
emission objects in common between adjacent FLAMES plate 
configurations. From this we determined the spectroscopic com¬ 
pleteness Cspec.fb- The success rates are affected by the spectro- 


^ https://www.eso.org/sci/software/pipelines/giraffe/giraf-pipe- 
recipes.html 

IRAF is distributed by the National Optical Astronomy 
Observatory, which is operated by the Association of Universities 
for Research in Astronomy (AURA) under cooperative agreement with 
the National Science Foundation. 

^ In IStrader et all ll201lh the ’nominal’ success rate is defined as the 
ratio between number of slit/fibres allocated to candidates and the num¬ 
ber of identified GCs. 


scopic completeness as well as by the presence of stars in the cat¬ 
alogue of candidates to which the fibres are allocated. In Table|2] 
we report the total number of allocated fibres together with the 
nominal success rate, the number of spectroscopically confirmed 
PNs, and the spectroscopic completeness Cspec.tb for all of the 
FLAMES fields. 

Because success rates are also dependent on the absence 
of stars in the submitted catalogue, candidates from the more 
stringent automatic sample have higher success rates than candi¬ 
dates from the visual catalogue, reaching > 70% in M87SUB1 
BrightFOS and M87SUB2 BrightFOl (see Table|2]i. On the other 
hand, about 30% of the confirmed spectra come from targets in 
the visual c atalogue: this is consisten t with the results of the sim¬ 
ulations by iLongobardi et alJ (1201 3h for the fraction of missed 
true line emission sources with automatic selection criteria. 

2.4. Classification of the extracted spectra 

The colour selection criteria are based on the strong [OIII] 
/15007A emission of a PN, with faint or no continuum. 
Nonetheless, background galaxies like Lya emitters at z ~ 3.1 
and [Oil] /13727.26A emitters at z ~ 0.34 have relatively strong 
lines that fall within the bandpass of the narrowband filter. Thus, 
we classify the extracted spectra on the basis of the shape of the 
line profile of the strongest emission. The extracted spectra fall 
into these categories: 

- PN spectra: the [OIII]/15007A emission of a PN is charac¬ 
terised by a narrow and symmetric line shape and very low 
continuum. In high S/N spectra, we detected the redshifted 
[OIII] /14959/5007A doublet. Typical S/N ratios for the spec¬ 
troscopically confirmed PN [OIII]]/15007A cover a range of 
2.5 < S/N < 8.5 per resolution element. In Fig. [^we show 
examples of single PN spectra with different Vlos- 

- Lya spectra; the emission line of a Lya emitter has a broader 
and more asymmetric line profile, characterised by a steep 
drop-off at bluer wavelengths. This kind of signature comes 
from the forest absorption bluewards of Lya: the symmetric 
emission line is truncated below the object redshift by Lya 
scattering in the intergalactic medium. Fig. [3] shows an ex¬ 
ample of an extracted FLAMES spectrum for a Lya emitter. 

- [Oil] spectra: the [OII]/l3727A emitters are characterised by 
the redshifted, resolved, and broad emission lines of the oxy- 
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Table 2. Summary of the total number of allocated fibres, nominal success rate, detected number of PNs, and spectroscopic completeness in all 
fields. Numbers in brackets refer to sources in common with the automatic sample. 


FLAMES 

Conf. 

# of targets with 
fibres allocated 

Nominal success rate 

Confirmed PNs 

Cspec,fb 

M87SUB1 Bright FOl 

33 (12) 

52% (42%) 

8(3) 

0.6 

M87SUB1 Bright F02 

41 (19) 

23% (32%) 

4(4) 

0.5 

M87SUB1 Bright F03 

33 (14) 

42% (71%) 

8 (8) 

0.9 

M87SUB1 Norm FOl 

125 (58) 

47% (57%) 

49 (26) 

0.9 

M87SUB1 NormF02 

130 (65) 

52% (62%) 

55 (35) 

0.9 

M87SUB1 NormF03 

127 (56) 

53% (64%) 

54 (31) 

0.9 

M87SUB2 Bright FOl 

23 (8) 

61% (75%) 

5(3) 

0.8 

M87SUB2 Bright F02 

31 (16) 

26% (25%) 

5(4) 

0.6 

M87SUB2 Bright F03 

26 (7) 

54% (57%) 

5(3) 

0.8 

M87SUB2 Norm FOl 

104 (40) 

16% (17%) 

10 (6) 

0.3 

M87SUB2 Norm F02‘' 

144(71) 

42% (51%) 

44 (27) 

0.6 

M87SUB2 Norm F03 

117(50) 

24% (30%) 

23 (14) 

0.4 

M87SUB CEN 

130 (94) 

54% (60%) 

63 (51) 

0.9 

M87SUB FEDGE 

131 (60) 

27% (30%) 

30(16) 

0.8 


Notes. Fibre configuration was modified between different exposures. 


gen doublet at /13726-3729A. Fig. |3] shows an example of 
an extracted FLAMES spectrum for an [Oil] line emitting 
galaxy. 




4960 4980 5000 5020 5040 

Wavelength (Angstrom) 


Fig. 3. Spectra for background emission line galaxies: Lya emitter at 
z ~ 3.1 (top panel) and [Oil] emitter at z ~ 0.34 (bottom panel). 
Red vertical lines mark the positions of the two components of the 
/13726+3729A emission. We smoothed the spectra to 0.015 nm per 
pixel. 


Fig. 2. Spectra of single confirmed PNs. The top panel shows the spec¬ 
trum for a PN dynamically hound to the halo component. The middle 
and bottom panels show spectra for PNs dynamically unbound to the 
halo (see Sect|3), with higher (middle panel) and lower (bottom panel) 
velocity than the M87 systemic velocity. Red vertical lines mark the 
positions of the two oxygen lines at their redshifted wavelengths. We 
smoothed the spectra to 0.015 nm per pixel. 


same fraction of the automatic catalogue is in fact La or [Oil] 
emitters. 


The final sample of emission line objects consists of 380 
sources, of which 287 were classified as PNs and the remain¬ 
ing as background emission line galaxies, either as Lya or as 
[Oil] emitters. This is the largest sample of spectroscopically 
confirmed PNs around M87 thus far, which is about a factor 15 
larger than the sample of iDohertv et al.l (l2009l) . 

The fraction of background emitters is con sistent within one 
cr with the estimate in the photometric study of iLongobardi et all 
(l2013h . They estimated that ~ 25% of the total imaging sample 
would be background emission galaxies. Here we find that the 


2.5. Accuracy of the velocity measurements 

From the repeat observations of the same candidates in ar¬ 
eas where different FLAMES plate configurations overlap (see 
Fig. [D, we obtained independent velocity measurements for 
a subsample of PNs. The median deviation of these measure¬ 
ments is 4.2 kms ' and the whole distribution covers a range of 
0.6 < AVlos ^ 16.2 kms '. The largest errors occur when acos¬ 
mic ray falls near the wavelength of the [OIII]/I5007A emission 
in one of the exposures. 
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3. Halo and 1C PN components 

When studying the outer regions of M87 out to a distance of 
~ 150 kpc from the galaxy centre, we are tracing the light 
in the radial range w here the M87 halo blends into the ICL. 
lArnaboldi et alJ (l2004l) showed that the M87 stellar halo and the 
Virgo core ICL coexist fo r distan ces > 60 kpc from the galaxy 
centre, and iDohertv et alJ (12009) showed that the two compo¬ 
nents overlap out to 150 kpc. iLongobardi et al.l (1201 3h showed 
that the observed slope of the PN number density profile is con¬ 
sistent with the superposition of two PN populations associated 
with the M87 halo and ICL, respectively. 

In the following subsections, we show that the distribution 
of LOS velocities obtained from the FLAMES spectra shows 
evidence for two dynamically distinct PN populations at large 
radii, conhrming this interpretation. 


3.1. Projected PN phase-space diagram 


R (kpc) 

0 50 100 150 200 



0 500 1000 1500 2000 2500 

R (orcsec) 


Fig. 4. Projected phase-space diagram showing Vlos vs. major 
axis distance R from the centre of M87, for all spectroscopi¬ 
cally confirmed PNs (black asterisks). The major axis distance 
is given both in arcsec (bottom axis), and in kpc (top axis), 
where 73 pc = 1". The dotted horizontal line shows M87’s sys¬ 
temic velocity Vsys = 1275 kms^' from Sect l3.2l The shaded 
area represents the region of the projected phase-space where 
the blue-shifted [OIII] 4959A emission line would fall below the 
wavelength of the blue edge of the FLAMES sort ordering hlter 
HR08. 

For each conhrmed PN spectrum, we measured Vlos and 
computed the major axis distance via the formula - 
xpN^/(l - e)^ -H Vn^T, wh ere e is the isophote’s ellipticity from 
iKormendv et al.l (l2009l) and XpN, ypN are the PN coordinates 
measured in a reference frame centred on M87, where the y axis 
is aligned with t he major axis of the o uter elliptical isophotes 
at P.A. = -25.6 (IKormendv et al.ll2009l) . In Fig.|4] we show the 
projected phase-space diagram Vlos.pn vs. R for the spectro¬ 
scopically confirmed PNs in the M87 survey (black asterisks). 

In this projected phase-space, PN velocities show a concen¬ 
tration around the systemic velocity of M87( Vsys = 1275kms“*; 
see Sect. 13.21) . in addition to a scattered distribution at higher and 
lower velocities. In Fig.|5j we show the histogram of the veloci¬ 
ties for the entire sample. It has two very strong and quite asym¬ 



-3000 -2000 -1000 0 1000 2000 3000 

Velocity (Km/s) 

Fig. 5. Histogram of the line-of-sight velocities of the spectroscopically 
confirmed PNs (black histogram) fitted with a double Gaussian (black 
curve). Red and blue lines represent the two Gaussians associated with 
the M87 halo and the IC components. 


metric wings around the main peak. The wing (or tail) at low 
velocities is more extended than the tail at high velocities with 
respect to the peak. These extended tails are very different from 
those measured in the LOS velocity distributions (LOSVDs) of 
isolated early-type galaxies, which are well de scribed by sin¬ 
gle Gaussian distri butions to within ~ 1% (e.g.. lGerhardlll993l: 
iBender et al.lll994h . The total velocity distribution is well fitted 
by a double Gaussian, while a single Gaussian is a poor fit to 
the observed LOSVD. The reduced of a double Gaussian fit 
is 1.1, while for a single Gaussian it is 2.0. Therefore, we fit the 
LOSVD in Fig. 5 with the sum of two Gaussians. The ht to the 
data returns a narrow component centred on Vlos.h = 1270.4 
kms^' with velocity dispersion of cr„ = 298.4 kms^*, and a 
broad component, centred on VLos.b = 999.5 kms * with a larger 
velocity dispersion cTb = 881.0 kms“'. 

The broad component is shifted from the M87 systemic LOS 
velocity; both Vlos, b and cTb are consistent with those values de¬ 
termined for the LOSVD o f galaxies in the main sub-cluster re¬ 
gion A of the Virgo cluster (iBinggeli et al.lflO^IConselice et al.l 
1200 ll) . The LOSVD around M87 is thus bimodal, containing a 
narrow component associated with the systemic velocity of the 
galaxy plus a broader component associated with the ICL. 

Different LOSVDs for the halo and ICL are predicted by 
cosmological analysis of structur e formation. Using hydrody- 
namic al cosmological simulations. lDolag et al.l(l^lol) : ICui et al.l 
(12014t) study the LOSVD of star particles at the centre of their 
clusters. These authors hnd that this LOSVD is well described 
by the sum of two Gaussians, with similar average velocities 
but different cr. One component is gravitationally bound to the 
galaxy and more spatially concentrated; the other is more diffuse 
and its high-velocity dispersion reflects the satellites’ orbital dis¬ 
tribution in the cluster gravitational potential. It is plausible that 
the halo stars are spatially confined as a conseque nce of merging 
processes that led to the formation of the BCG (iMurante et al.l 
120071: IContinreF^l20l4 . 

3.2. Robust separation of the M87 haio and ICL 

The overall LOSVD of the PNs in M87 is characterised by a 
narrow component associated with the M87 halo, superposed on 
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a broad IC component with a shifted mean velocity and much 
larger velocity width. We are now interested in separating halo 
and IC PNs based on their different L OS VPs. In this a nalysi s, we 
combine our PN sample with that of iDohertv et al.l (l2009ll . We 
concentrate on those PNs that have a major axis distance R < 
190 kpc to study the transition between the M87 halo and the 
Virgo ICL. The combined total sample consists of 299 PNs with 
meas ured positions and ve locities within R < 190 kpc; the PNs 
from iDohertv et al.l (l2009l) further out are classified as ICPNs 
and are not discussed further. 

To separate halo and IC components, we use a sigma clipping 
algorithm in elliptical radial bins, corresponding to vertical strips 
in the projected phase-space diagram in Fig. |6] The idea is to 
separate the velocities in the narrower Gaussian (see Fig.|5]l from 
those in the high- and low-velocity wings of the distribution. By 
tagging our sources depending on whether their Vlos belongs 
to the narrower or wider Gaussian, we assign PNs to either the 
M87 halo or the ICL. 

The velocity dispersion profile of the M87 halo: a robust 
sigma estimate- We binned the PN velocity sample in elliptical 
annuli and, for each bin, we determined the standard deviation 
of the LOSVD for the PNs in this bin. We applied a 2cr limit 
with respect to the systemic velocity of M8'0, Vsys = 1275 + 24 
kms ', and calculated the dispersion for all PNs with |Vlos- 
Vsysl < 2crkms”*. We then scaled this dispersion by a numerical 
factor determined from Monte Carlo simulations to correct to the 
dispe rsion of a complete Gaussian distribution (see lMcNeil et al.l 
l201(]li . Since we expect the initial estimate for the 2cr to be influ¬ 
enced by the ICPNs, we repeat this process until the dispersion 
value stabilises. 

Separating M87 halo and IC PNs: sigma clipping - We now 
identify the two components using a sigma clipping algorithm. 
We begin by classifying as ICPNs all velocity outliers that devi¬ 
ate from the M87 systemic velocity Vjys = 1275 kms^' by more 
than 2cr. To obtain the required cr value at the radius of each 
PN, we took the robust estimates of the velocity dispersions in 
the elliptical radial bins, and fitted these data with a fourth-order 
polynomial. This takes radial gradients into account and, at the 
same time, reduces the effects of binning and scatter in the dis¬ 
persion profile on the separation of the components. Using the 
2cr threshold from the interpolated polynomial, we identify 243 
objects as M87 halo PNs. 

Two further steps are still needed. Firstly, the 2cr criterion 
accounts for ~ 95.5% of a complete Gaussian distribution; hence 
we expect it to have missed 11 halo PN candidates. To include 
these, we considered all outliers within 3cr from the M87 Vsys 
and, from those, selected the 11 with the smallest |VLos-Vsysl/cr 
ratios. This leads to a final sample of 254 M87 halo PNs and 45 
ICPNs. In this final s ample, 11 halo PNs and 1ICPN came from 
IDohertv et al.l (l2009l) . 

Secondly, as can be seen in Fig. |5] the ICL and halo veloc¬ 
ity distributions overlap, and as result, when using the sigma¬ 
clipping algorithm, the ICPNs at low velocities relative to Vsys ~ 
1275 kms ' are considered to be part of the halo component. To 
statistically quantify this effect, we compared the halo and IC 
velocity distributions in each radial bin (or slice in the phase- 
space). We approximated the halo distribution in the bin as a 
Gaussian centred on the systemic velocity of M87 with a dis- 


® The systemic velocity of M87 is taken as the median value of 
the entire sample of velocities within two sigma of the original me- 
dian. Our valu e 1275 ± 24kms“' is consistent with t hose obtained by 
iBinggeli et all ( 119871) (Vsys = 1258 ± lOkms^') , and ICappellari et ^ 
(1201 Ih tV„,. = 1284 ± 5kms-‘). 


pension equal to the average value for that bin, and for the IC 
component we used the same Gaussian in all radial bins using 
the same parameters as in Fig. |5] We then calculated the frac¬ 
tion of ICPNs that lie inside the halo distribution as the area of 
overlap between the two curves. With this analysis we obtain a 
statistical estimate of the number of ICPNs contained in the halo 
velocity distribution in each bin. For all bins combined, we find 
that a further ~17% of the halo sample, i.e., ~ 44 PNs, are to be 
associated with the IC component. 

In Fig. |6] we show the projected phase-space distribution of 
our PNs as in FiglU with velocities colour-coded to show the 
membership to the halo (red) and ICL (blue). Since we only 
know statistically, but not individually, which PNs in the halo 
velocity range are ICPNs, these are also shown with the red halo 
colour. 

In Fig. |6]we also show 1/2/2.5cr limit contours for the halo 
PNs, obtained by fitting a polynomial to the cr values from the 
robust estimation. The velocity dispersion profile increases from 
250" to 1200 ", and then decreases, showing a colder component 
at radii R > 1200". A more detailed analysis of the dispersion 
profile will be provided in a separate paper (Longobardi et al., in 
preparation). 

We now enquire whether any ICPNs could be associated 
with other galaxies in the cluster. The first to consider is M86, 
because it is bright, has the most extended halo, and is rela¬ 
tively close_to_^487_(~ 1.2° a way). Extrapolating the Sersic fit 
from iJanowiecki et al.l (1201 Ol) . we compute the total luminos¬ 
ity associated with the M86 halo in our surveyed area to be 
Lv.boi ~ 8.0 X lO^Leboi- We adopt a similar PN specific num¬ 
ber as for M87, 02.5 = 1.01 x lO^^PN/Lo^oi (see Sect l5.1b . but 
we correct for the fact that we can only sample the brightest 
2.3 mag of the PNLF because of the larger distance to M86 
(iMei et al.ll2007h . This leads to a predicted number of 7 M86 
PNs that could be counted as IC component, i.e., only ~ 2.5% 
of the completeness-corrected ICPN sample (see Sect l5.11 l. or 2 
out of a total of 88 ICPNs. 

Other fainter galaxies have less extended halos and con¬ 
tribute less to the sampled ICPNs around M87 unless they are 
very near or in the survey area. Moreover, from the shape of 
the ICPNLF (see Sect l5.21 l the majority of ICPNs is likely to 
have been stripped from low-luminosity star-forming galaxies, 
such as M33 and the LMC. Any ICPNs still bound to dwarf 
galaxies would be correlated with their position and velocity. 
Only three PNs with velocity V~350 km/s out of the ICPN 
sample seem correlated with the position and systemic veloc- 
ity of a dwarf galaxy, IC3549, w hich has Vlos = 375 kms^' 
(lAdelman-McCarthv & et al. 20091) . Hence, we can state that the 
contribution from other nearby galaxies to the ICPN sample is 
negligibile and does not affect any of the results we present. 

Finally, a few ICPNs are characterised by extraordinary blue 
shifts relative to M87 (Vlos < -lOOOkms^*). These (hyper-) 
velocities of |VLos-Vsysl > 2300kms^* relative to M87, corre¬ 
sponding to several cryirgo, could be because of infall from the 
outskirts of the Virgo cluster, perhaps associated with the infall 
of the M86 group, or they could be tracers of three-body interac¬ 
tions, as previously hypothesised for the extreme globular cluster 
observed at a projected distance of ~ 80 kpc from M8 7 with a 
velocity of Vlos < -1025 kms ' dCaldwell et al.ll2()T4l) . 

3.3. Spectroscopic validation of the PN subsample 

In Section IZ4l we classified 287 spectra as PNs on the basis of 
the line profile of the strongest emission. We now strengthen this 
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Fig. 6. Projected phase-space diagram, Vlos vs. major axis dis¬ 
tance from the centre of M87, for all spectroscopic ally con¬ 
firmed PNs from this work and iDohertv et al.l (l2009l) .The ma¬ 
jor axis distance is given both in arcsec (bottom axis), and in 
kpc (top axis), where 73 pc = 1". The PNs are classified as M87 
halo PNs (red asterisks) and ICPNs (blue asterisks), respec¬ 
tively; see text. Black squares identif y spectroscopically con¬ 
firmed PNs from iDohertv et all (l2009l) . The smoothed 1, 2, and 
2.5 cr thresholds are represented by the dashed, dotted, and dot- 
dashed lines, respectively. The dashed horizontal line shows the 
M87 systemic velocity Vsys = 1275 kms“* as computed in 
Section lT2l while the continuous gr een line shows the effec tive 
radius = 703.914" determinedbv iKormendv et al.l (l2009l) . At 
Vlos = -220 kms^*, we plot the M86 systemic velocity (long 
dashed line). 

earlier classification based on the detection of the weaker 4959A 
line of the [OIII] doublet in the PN spectra. 

In 114 out of the 287 spectra, nearly 40% of the sample, we 
are able to detect the Doppler-shifted [OIII] 44959A line with 
the expected 1:3 ratio of the [OIII]/14959/5007A line fluxes. 
In the rest of the sample, the single spectra do not have the 
required S/N (S/N > 3) for the main [OIIIj/15007 A emission 
line to all ow the detection of the weaker 4959A line. However, 
following lArnaboldi et alJ (l2003h we can achieve the required 
S/N by stacking these spectra, after shifting their [OIII]/15007A 
emission to a common reference wavelength. By measuring the 
[OIII]/15007/4959A line ratio of the coadded spectrum, we can 
statistically constrain the fraction of misclassified PN spectra: if 
the stacked spectrum contains misclassified PN candidates, the 
ratio [OIII] /15007/4959A is larger than three. 

The 287 spectra are further grouped in three classes: M87 
halo spectra, IC high-velocity, and IC low-velocity spectra. We 
have 254, 13, and 31 PN in each category, respectively. Three 
sources in the IC low-velocity class are not included because 
their Doppler shifted [OIII]/l4959A emission fall at a shorter 
wavelength than the blue edge of the FLAMES HR08 filter (see 
Fig. m. The spectra are shifted so that the main [OIII]/15007A 
emission falls at the same nominal wavelength for all spectra of 
a given class. For the halo and IC high-velocity class, we adopt 
the nominal wavelength of 5029A for [OIII]/15007A, i.e. the red- 
shifted wavelength at the systemic velocity of M87. For the IC 
low-velocity class, the adopted nominal wavelength is 5000A. 
We shifted the single spectra for each class, and obtained three 


coadded spectra. These are plotted in Fig. |2l The [OIII] dou¬ 
blet is visible in all spectra, with the [OIII]/14959A line visible 
at the correspondent shifted wavelength. In all three spectra, the 
FWHM of [OIII]45007A is FWHM^O.hA (Table [3]l, somewhat 
smaller than the typical FWHM of single spectra, FWHM=0.8 
A. Because of filter edge effects, together with the shifting and 
alignment of the spectra with reference to the stronger line, the 
weaker [OIII]/l4959A emission is broader in the IC low-velocity 
spectrum than in the single PN spectra. 

The fluxes and FWHMs of the [OIII]44959/5007A emis¬ 
sions in the coadded spectra are measured via a Gaussian fit to 
the lines with the IRAF task splot . The errors on the line fluxes 
were calculated using the relation dPerez-Montero & Dfazll20()^ 

0-1 = + EW/{NA)]^^^, ( 1 ) 

where cri is the error in the line flux, cr^ is the standard deviation 
in a box near the measured line and represents the error in the 
definition of the continuum, N is the number of pixels used to 
measure the line flux, EW is the equivalent width of the line, 
and A is the wavelength dispersion in A pixel '. In Table [3] we 
give a summary of the observed properties of the [OIII] doublet 
of the coadded spectra. 

The [OIII] doublet is visible in the three spectra from the 
different classes, with approximate flux ratio of 3:1, as expected 
from atomic physics. Given the uncertainties on the measured 
fluxes, the spectral validation shows that at Icr > 99% of the 
halo and > 90% of the IC high-velocity spectra are true PNs. 
This implies that only two halo PNs and one high-velocity ICPN 
could be misclassified. For the IC low-velocity PNs, the line ratio 
of the [OIII]5007/4959 A is 3.1 + 0.3 from Tabled Hence up to 
13% of the spectra could be misclassified, i.e. four ICPNs. On 
the total sample of 287 spectroscopically confirmed PNs, at most 
seven candidates might have been misclassified, i.e. only 2% of 
the whole sample. 

4. Halo and IC density profiles 

By tagging PNs according to their Vlos, we differentiated be¬ 
tween halo and IC PNs in the last section. We are now interested 
in recovering the spatial distribution of these two components 
and in studying their number density profiles separately. 

In Fig. [8] we show the sky positions of the spectroscopically 
confirmed halo and IC PNs. If the halo and ICL are characterised 
by different evolutionary histories, their radial distributions may 
be different. If so, this can be seen in the PN density profile. 
Because PNs follow light, the presence of a single parent stel¬ 
lar population would be reflected in a PN number density dis¬ 
tribution proportional to the surface brightness profile at each 
radius. Deviations of the PN number density from the light pro¬ 
file would trace the pre sence of more than one stellar population 
in the surveyed region (iLongobardi et al.ll2013l) . 

4.1. Completeness corrections 

To construct a PN density profile, we require a spatially complete 
sample hence we must correct the detected PN number using a 
completeness function, Ctot(xpN, ypN), which accounts for the se¬ 
lection function of the sampled PNs over the surveyed area. Our 
PN sample is affected by four different kinds of incompleteness, 
which are related to the photometric identification of the candi¬ 
dates and the selection effects in the spectroscopic observations. 
These are: 
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Table 3. Fluxes, FWFlMs, and line ratios of the [OIII] /15007/4959A doublet for halo, high-velocity and low-velocity IC PNs sub-samples. 


IDs PN sub-sample 

Line FIUX 5007 
( Counts ) 

( 1 ) 

Line FIUX 4959 
(Counts ) 

( 2 ) 

FWHM 5007 

(A) 

(3) 

FWHM 4959 

(A) 

(4) 

Line Ratio 

(5) 

halo 

1528.0±0.2 

506.9±0.3 

0.60±0.08 

0.76±0.1 

3.01±0.02 

IC high-velocity 

64.1±3.0 

20.6±L2 

0.60±0.07 

0.6±0.09 

3.1±0.2 

IC low-velocity 

106.3±3.0 

33.9±3.7 

0.60±0.08 

L4±0.1 

3.1±0.3 


Notes. Columns (1) and (2): fluxes for the [OIII] main and second emission line, respectively. Errors are calculated following Eq[T] 

Columns (3) and (4): emission line FWHMs from Gaussian htting of the spectral lines. Errors represent Icr uncertainties from sampling statistics. 
Column(5): flux ratio between the [OIII] doublet emission lines. Errors are propagated from the flux errors. 


Photometric incompleteness, characterised by; 

i) Cphot,sp(^) - Spatial incompleteness due to the high 
galaxy background and foreground stars that affect the 
detectio n of PN candidates in the images. This was esti¬ 
mated in lLongobardi et akl (1201 3h by adding a simulated 
PN population to the scientific images and determining 
the fraction of simulated objects recovered by SExtractor 
in the different elliptical annuli shown in Fig. [8] 

ii) Cphot,coi - Colour incompleteness due to the colour crite¬ 
ria adopted for the automatic selection of PN candidates. 
This colour incompleteness was computed by analysing 
the properties of the recovered simulat ed PN population 
in the colour-magnitude diagram (see iLongobardi et al.l 
1201 3l for more detail). This incompleteness affects only 
the M87SUB FCEN and M87SUB FEDGE fields, where 
the photometric candidates were selected through colour 
criteria (see Sect. IHl, and is magnitude-dependent. For 
constructing the density profile, however, we use an av¬ 
erage value, computed for the whole sample down to 2.5 
magnitudes below the bright cut-off, which amounts to 
0.7. 

Spectroscopic incompleteness, characterised by: 

iii) Cspec,sp(f^) - Spatial incompleteness due to the limited 
number of fibres (up to 132) that can be allocated for 
each FLAMES field. This was estimated by computing 
the ratio between the number of allocated fibres and the 
total number of photometric candidates in each elliptical 
annulus shown in Fig.|8] 

iv) Cspec,fb(xpN, Ypn) ' Incompleteness due to fiber-target 
misalignment. This incompleteness was estimated from 
the detection statistics of objects in common between 
overlapping FLAMES plate configurations, whose [OIII] 
emission was detected in spectra taken with either one or 
two of the two-plate configurations (see Sect. I23 T i. 

Our spectra are sufficiently deep that we found no depen¬ 
dence of the spectroscopic incompleteness on the magnitude 

of the PN candidates. 

The total completeness function, Ctot(xpN, ypw). is the prod¬ 
uct of the photometric and spectroscopic incompleteness: 

Otot(XpN5 Ypn) ~ Cphot,sp(^) * Cphot.col 

* Ospec,sp(f^) * Cspec,fb(XpN, Ypn), (2) 

with 

_ / Cphot.coi for M87SUB FCEN, M87SUB FEDGE, 
phot,col “ i 1 elsewhere. 


The completeness-corrected number of PNs in each bin of major 
axis distance R is then 


N,(R) =Xf*k,(R), (3) 

where the sum extends over all PNs of the halo or IC compo¬ 
nent in the bin, respectively, and k,(R) = 1 /Ctot(xpNi, ypNi) is the 
completeness-corrected specific weight of each observed PN at 
its position. 


4.2. Density profiles of halo and IC component 

To construct the PN number density profile and compare it to 
the galaxy’s surface brightness profile, we bin our PN sample in 
elliptical annuli. The radial range of the elliptical annuli is cho¬ 
sen such that they include the major axis distance of the inner- 
most and outermost PN candidates in the photometric sample of 
ILongoba rdi et al.l (l2013h. Th eir PA.s and ellipticities are taken 
from iKormendv et al.l (1200^ . The sizes of the annuli are deter¬ 
mined separately for the halo and IC components, such that for 
each component all bins contain at least ten spectroscopically 
confirmed PNs. In each annulus and for each component, we 
compute the completeness-corrected PN number density as the 
ratio of the completeness-corrected number of PN (Eq. [3 and 
the area of the portion of the annulus intersecting our FOVO, 


O'PNj(f^) - 


A(R) ■ 


(4) 


Here we consider all the spectroscopically confirmed PNs whose 
magnitudes are within 2.5 mag below the bright cut-off. The 
subscript j indicates the two different PN components, halo and 
ICE. The major axis distances R in Eq.|4]are determined by com¬ 
puting the average major axis distance of all PNs falling within 
each elliptical annulus. 


In Fig. 0we show the comparison between the V-band M8 7 
surface brightness profile pkoq from iKormendv et al.l (l2009ll . 
with the logarithmic PN density profiles for the halo (red tri¬ 
angles) and IC (blue diamonds) PNs, defined as 

PPNj ( R ) = -2.5 logio (crpN.(R)) -H po - (5) 

The value po is a constant to be added so that the PN number 
density profile matches the pko 9 surface brightness profile. As 
described in Sect. 13.21 the kinematic decomposition of the halo 
and IC components does not identify ICPNs in the velocity range 

^ These areas are estimated using Monte Carlo integration. 
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Fig. 7. Combined spectra for the spectroscopically confirmed sub¬ 
samples of PNs belonging to the halo (upper panel), IC high-velocity 
(central panel) and IC low-velocity (bottom panel) classes. As expected 
from atomic physics, the [OIII]/i4959/5007A doublet is detected with 
a flux ratio 1:3 confirming that the large majority of our candidates are 
true PNs (see text for more details). 


of the halo, and their contribution must be evaluated statistically 
for each bin. Averaged over the bins, ~10% of the halo PN sam¬ 
ple is thus estimated to be associated with the ICL, in addition 
to the ICPNs identified from their large velocities. For each of 
the radial bins, we subtract the estimated contribution from the 
halo component and add it to the IC component. In Fig. |9] the 
PN number densities for halo and ICL account for this effect. 
Furthermore, the profiles shown in Fig. 0 are computed using 
all PNs whose magnitudes are within 2.5 magnitude from the 
bright cut-off. However, the halo and IC number density profiles 
do not change significantly when the fainter spectroscopically 
confirmed PN are also included. 

We find that the agrees well with the surface photom¬ 

etry: the halo PN logarithmic number density profile follows the 
surface brightness profile. The ICPN logarithmic number den¬ 
sity, as the halo profile, is centrally concentrated towards M87. 
However, it has a flatter profile that decreases towards larger 
radii as /icl with y in the range [-0.34; -0.04], depend- 



188.2 188.0 187.8 187.6 187.4 187.2 

RA (degree) 


Fig. 8. Sky positions of the spectroscopically confirmed halo 
(red crosses) and IC PNs (blue crosses). The green plus sign indi¬ 
cates the centre of M87. The ellipses (dotted lines) trace the M87 
isophotes between photometric major axis ra dii R - 2.'8 and 
R - 40 .'7, for position angle PA.=-25.6°, from iKormendv et al.l 
(l2009h . The solid squ ares depict the area cove red by our narrow- 
band imaging survey (iLongobardi et al.ll2013h . North is up, East 
to the left. 



Fig. 9. Logarithmic PN number density profiles for the M87 halo 
(red triangles) and the IC PNs (blue diamonds), corrected for 
incompleteness. The error bars include counting statistics and 
completeness correction. The halo PN profile compares well 
with the surface brightness profile from iKormendv et al.l (l2009l 
plus signs). The continuous black line shows their Sersic fit 
with n = 11.8. The ICPN number density decreases towards 
larger radii as Ijcl R^ with y in the range [-0.34; -0.04] 
(full and dashed-dotted blue lines). Black squares show the com¬ 
bined halo and IC PN number density profile, which is well 
modelled by the two-component photometric model (dot-dashed 
black line, see Section|5]l. 
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Fig. 10. Top panel: Same as in Fig. |3 but here the halo and IC 
number density profiles are computed separately for the NW 
(filled symbols) and SE (open symbols) sides of M87. Bottom 
panel: Difference between the logarithmic PN number density 
profiles of the halo (red dots) and IC (blue dots) PNs in the NW 
and SE M87 regions. A value of zero indicates no difference in 
the number of sources on opposite sides at a given major axis 
distance. 


ing on the choice of binning. These results are consistent with 
predictions from hydro-dynamical simulations, where the radial 
density profile of the bound component, i.e. the halo, is ob- 
served to be much steeper than that of t he diffuse IC component 
dMurante et al .11200^ iDolag et^l2010l) . 

Einally, we compute density profiles for the PNs for the NW 
and SE sides of M87 independently, to search for a NW-SE 
asymmetry in the spatial distribution. Eor this test, we now use 
the total sample of PNs, including very faint PNs, to increase 
the statistics in each half annulus (this also changes the bin¬ 
ning). Pig. [To] shows the four number density profiles; each pair 
of profiles are consistent with each other and the halo profiles 
are consistent with the galaxy surface brightness profile. Thus 
within the statistical uncertainties, the stellar halo density is NW- 
SE symmetric and the PN number density follows the light on 
both sides. The ICPN density profile pair are flatter than the halo 
again, and consistent with the previous surface brightness profile 
IiCL with y - [-0.34; -0.04]; no asymmetry is evident. In 
Pig.dOjwe also plot the differences between the NW and SE log¬ 
arithmic densities as function of radius (bottom panel). Eor both 
the halo (red filled circles) and the ICE (blue filled circles), these 
differences are consistent with zero within the uncertainties (see 
further discussion in Sect.lbj). 

5. Halo and IC populations 

5. 1. The a parameter 

The total number of PNs, Npn, is proportional to the total bolo- 
metric luminosity of the parent stellar population through the 
luminosity-specific PN density, or a-parameter for short, such 
that Npn = crLboi- The a-parameter determined for this PN sam¬ 
ple is an estimate of the total number of PNs within 2.5 mag¬ 
nitudes of the bright cut-off of the PNLE, because of the mag¬ 
nitude limit of this survey at m5007 =28.8 and the bright cut-off 


at m*^Q^y^-2()3 (see Section fS^ . The measured value of a is de¬ 
rived from the scaling factor required to match the PN number 
density profile to the surface brightness profile in the V band and 
then taking the appropriate bolometric correction into account. 

In Pig. |9| we showed the total PN number density profile 
for the spectroscopically confirmed PNs sample within 2.5 mag¬ 
nitudes from the bright cut-off, together with the halo and IC 
PN number densities. The total PN number density profile flat¬ 
tens at large radii compared to the V-band su r face b rightness 
profile. As already shown in iLongobardi et alJ (l2013h . this ef¬ 
fect is due to a superposition of stellar populations with different 
PN specific frequencies. The flattening for the spectroscopocally 
confirmed PN number density pro file provides independent sup¬ 
port for these results. Poliowing ILongobardi et alJ d2013h . we 
can model it by two components such that 


d-(R) - [ahaloI(f^)halo,bol + Q'ICLl(f^)lCL,bol] 


“ ^^halo 


I(f^)K09,bol + 


/OICL 

V^halo 


- 1 


I(^)lCL,bol , 


(6) 

(7) 


where d-{R) represents the predicted total PN surface density in 
units of NpNpc^^; I(/?)haio and I(f?)icL, with and without sub¬ 
script bol, are the bolometric and V-band surface brightnesses 
for the halo and the IC components, in Lopc^^, respectively; Iko 9 
is the M87 lum i nosity profile in the V band in LoPC^^, from 
iKormendv et al.l (l2009l) . The surface brightnesses I(/?)baio and 
I(f^)iCL are given by the Se rsic fit to the observed M87 surface 
brightness data (n = 11.8. IKormendv et al.ll2009l) . and by the 
scaled power-law fit to the IC surface density data from Sect. |4| 
respectively. They satisfy the relation Iko 9 = I(f^)baio + I(f^)iCL. 
which determines the normalisation of the IC surface brightness 
profile. In the surveyed area over a radial range 7kpc < R < 
150kpc, these profiles give total V-band luminosities Lhaio = 
4.41 X IO'^Lq and Licl = 0.53 x IO^^Lq, and after the bolo¬ 
metric correction (see below), total bolometric luminosities for 
the sampled halo and ICE of Lhaio,bo; = 9.05 x 10'°Lo,boi and 
Licl, bo; = 1.1 x 10^°Lo_boi- 

The surface luminosity d-{R) can be related to the bolometric 
surface brightness through the formula 


p{R) = -2.5 logio d-{R) + po, (8) 

where po is given by the analytical function 


Po - 2.5 logjg ahaio + K -I- (BCq - BCv). (9) 

In Eq.|2 tthaio is the specific PN number for the halo, K = 26.4 
mag arcsec'^ is the V-band conversion factor from mag arcsec”^ 
to physical units Lopc^^, BCq = -0.07 is the V-band bolometric 
correction for the Sun, and BCv=-Q35 is the bolometric cor¬ 
rection for the V-band (iBuzzoni et al.ll206^ . According to their 
simple stellar population (SSP) models for irregular, late, and 
early-type galaxies, this value can be used with 10% accuracy. 

Prom the offset value po - 16.56 + 0.08 mag arcsec^^ de¬ 
termined from the density profile in Pig.|9l we compute ahaio = 
(1.06 + 0.12) X 10-^ PN Prom Eqs.|6jand|2l the derived 
value for aicL is then (2.72 + 0.72) x 10“* PN using the 
steeper slope -0.34 for the ICE, but the difference for the shal¬ 
lower slope is only 5%. The difference in a-parameters is then 
~ 2.3cr. 

These lumino sity specific PN a value s are consistent with 
those obtained by ILongobardi et alJ (l2013h from the photomet¬ 
ric sample, and are now independently validated on the ba¬ 
sis of the spectroscopically confirmed PNs. Moreover, they 


11 





























A Longobardi et al.: Distinct stellar halo and intracluster light in the giant Virgo galaxy M87 


also agree withi n the uncerta i nties with the preoviously d eter¬ 
mined values bviDurrell et ^ (l2002h : iDohertv et alJ (l2009l) : see 
iLongobardi et alJ (I2OI3I) . The consistency of the spectroscopic 
and photometric values conhrms the accuracy of the estimated 
contamination by Lya background objects in the photometric 
sample. 

We can now compare our ahaio and cticl values with the 
known a values for PN populations in nearby galaxies. Galaxies 
with integrated {B - V) colours smaller than 0.8 are empiri¬ 
cally characterised by similar values of the a parameter, a ~ 
3 X 10“^NpnL“'|^^j, with a scatter of a factor of two. For red¬ 
der galaxies with {B - V) > 0.8, the spread of the measured 
values increases, spanning a range from a ~ 10“®NpN/Lo_boi to 
~ 6 X lO^^NpN/Lo boi- For these redder galaxies, there is an em¬ 
pirical inverse correlation of the a values with the FUV-V inte¬ 
grated colours of the parent stellar population, such that smaller 
a values are associated with galaxies with larger FUV-V excess. 
Hence, observationally, the a values are linked to the metallic- 
ity and st ar formation his t ory of the paren t stellar population of 
the PNs dPeimbertj 19901: iHui et al.ll l993t iBuzzoni et al.ll200^ 
ILongobardi et aklbool) . 

Our result that the IC component contributes more PNs per 
unit bolometric luminosity than the M87 halo light therefore sig¬ 
nals a change in the stellar populations from halo to ICL, con¬ 
sistent with the existence of a grad ient towards bluer colours of 
the M87 stellar light at large radii (iLiu et al.ll2005l: iRudick et alJ 
I 2 OIOI) . We interpret this gradient as the result of the gradual tran¬ 
sition from the redder halo light to the bluer ICL with decreasing 
surface brightness. 


5.2. The halo and ICL PNLFs 

The PNLF of [OIII]/15007A emission line fluxes is often 
empirically described via the truncated exponential formula 
dCiardullo et al.ll 19891) 

N(M) = {1 - , (10) 

where ci is a normalisation constant, C 2 - 0.307, and 

M*(5007)=-4.51 mag is the absolute magnitude of the PNLF 
bright cut-off. This analytical formula is designed; (i) to repro¬ 
duce the high-mass cut-off observed for the PN central stars 
in nearby galaxies, and ii) to model PNs as uniformly expand- 
ing homogeneous spheres io nised by a non-evolving central star 
(iHenize & Westerlundlll96^ . 

Our deep and extended imaging survey of PNs in the outer 
regions of M87 shows that the PNLF for this galaxy has sig- 
nihcant deviations from Eq. [TO] in the faint magnitude bins 
(ILongobardi et al.ll2013h : its slope ~l-2 magnit udes below the 
bright cut-off is steeper than expected from the ICiardullo et alJ 
(Il989t) formula. This is also true for the spectroscopically con- 
hrmed PN sample, as we show in Fig.[TT]in which we compare 
the total PNLF with the analytical formula for a distance modu¬ 
lus of 30.8. 

Our data allow us to analyse separately and compare the 
PNLFs of the two spectroscopically conhrmed PN samples for 
the M87 halo and ICL, from the bright cut-off down to 2.5 mag 
below. In the upper panel of Fig. [12] we show the PNLF for the 
spectroscopically conhrmed sample of halo PNs, corrected for 
detection incompleteness as a function of magnitude. The data 
points trace a smooth function, and the ht of the generalised an¬ 
alytical formula (Eq. ITOl) to the observed PNLE within the 2.5 
mag limit results in C 2 - 0.72 and a bright cut-off at i'n*^QQi = 26.3 


1000 


CD 

D 

E 


.0 100 

o 


10 


Fig. 11. Luminosity function of all spectroscopically conhrmed 
PNs corrected for i ncompleteness (full dot s). The dashed green 
line represents the ICiardullo et al.l (Il989l) formula for a dis¬ 
tance modulus of 30.8, convolved with photometric errors and 
normalised to the data at the bright end. The spectroscopi¬ 
cally conhrmed PNLE shows an excess of fainter PNs with re- 
spect to the analytical fo rmula, similar to the photometric PNLE 
(ILongobardi et alJl2013h . Open circles represent magnitude bins 
where the sample is not completeness-corrected. The error bars 
show the uncertainties from counting statistics and completeness 
correction. 
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(overplotted on the data). With M* = -4.51, this corresponds to 
a distance modulus (m - M)=30.8. 

In the central panel of Eig. [T2| we present the ICPNLE (full 
blue circles), corrected for incompleteness. As for the M87 halo 
PNLE, the ICPNLE is consistent with the same distance modulus 
as for M87. However, it shows an overall shallower gradient than 
the M87 halo PNLE: the best ht of Eq.[T0|to the empirical PNLE 
returns C 2 - 0.66. 

In addition to the shallower gradient at fainter magnitudes, 
the ICPNLE shows a clear dip at 1-1.5 mag fainter than 
This feature is statistically signihcant; the difference between the 
number of PNs in these bins with respect to the magnitude bins 
before and after the dip is > 3 o-^ip on both sides, where cr^ip 
is the uncertainty from Poisson statistics in the magnitude bins 
where the dip occurs. 

Dips in the PNLE are observed for PN populations de¬ 
tected in star-forming galaxies (irregulars/disks), and are ab¬ 
sent in the PNLEs of bulges or early-type galaxies. The mag¬ 
nitude below the bright cut-off at which the dip occurs varies 
between dif ferent PN populations, from ~ 2 to ~ 4 ma ; 
below M* (Ijacob v & De Marco 20021: Ciardullo et al. 200' 


iHernandez-Martmez & Pefial l2009t [Reid & Parker 2010 . We 
discuss this issue further in Section 16.11 In the bottom panel 
of Pig. [T2| we compare the PNLEs for the Virgo ICPNs and 
for the spectr oscopically conhrmed sample of PNs for the outer 
disk of M33 (ICiardullo et al.ll200^ . Both LEs are corrected for 
foreground Galact ic extinction, adopti ng reddening values of 
£(B-y)TrT,=0 02 (Ciardullo et al.lll998l) . and£'(B-y)M33=0.04 
(ICiardullo et al.l2004l) for the Virgo ICL and M33, respectivel}0. 
Absolute magnitudes are determined usi ng a distance modulu s 
of 30.8 for the ICL and 24.86 for M33 (ICiardullo et alJl200l . 
Both PNLEs show dips relative to the smooth luminosity func¬ 
tion: the ICL at 1-1.5 mag fainter than M*, the outer disk of M33 
at ~ 2.5 mag fainter than the bright cut-off. 

Pinally, we recall that about 17% of the PNs contained in the 
M87 halo PNLE shown in Pig. [T2|are ICPNs whose velocities 
fall in the same velocity range of the M87 halo PNs, and can 


® We used the relation of ICardelli et all dl989h with Ry = 3.1 to go 
from reddening to extinction at 5007A. 
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Fig. 12. Top panel: Luminosity function of the spectroscopically 
confirmed halo PNs, corrected for incompleteness (red circles). 
The red line shows the ht of the generalised analytical formula 
(Eq. [Tol l to the halo PNLF, with C 2 - 0.72 and bright cut-off at 
26.3 mag, corresponding to a distance modulus m - M - 30.8. 
Central panel: Completeness-corrected PNLF for the ICPNs 
(blue circles). The blue line shows the fit of the generalised an¬ 
alytical formula (Eq. [TOll to the ICPNLF, with C 2 - 0.66. The 
ICPNLF shows a dip 1-1.5 mag fainter than the bright cut-off. 
Bottom panel: PNLFs in absolute magnitude for the spectroscop¬ 
ically confir med Virgo ICPNs (blu e circles) and for the outer 
disk of M33 (ICiardullo et al.ll^04 cyan triangles) matched at 
the first bin. For this comparison, the data are corrected for line- 
of-sight reddening (see text). Absolute magnitudes are deter¬ 
mined us ing a distance modul us of 30.8 for the ICL and 24.86 
forM33 (ICiardullo et al.ll200l . 


therefore not be individually identified. The hint of a slight dip 
in the halo PNLF at ~ 1 mag below the bright cut-off may be 
due to these ICPNs. 

To summarise, the observed properties of the M87 halo and 
Virgo ICPN populations i.e. their a-parameters and PNLFs, 
show signihcant differences. Because these quantities depend on 
the physical properties of the parent stellar populations, these 
differences imply that the M87 halo and ICL consist of differ¬ 
ent populations of stars. To understand this better, more work 
is clearly required for a better theoretical understanding of how 
metallicity, age, and different star formation histories affect the 
post-AGB phases of stellar evolution and the resulting PN pop¬ 
ulations. 


6. Discussion 

6.1. The distinct haio and tC popuiations around M87 

In Sect. |3] we presented the projected phase-space distribution of 
the spectroscopically conhrmed PNs in our M87 fields (Fig.|6l). 
With a robust procedure we showed that the PN velocity distribu¬ 
tion splits into two kinematically very different components: the 
M87 halo (with mean velocity VLos.n = 1275 kms^' and veloc¬ 
ity dispersion (Xn - 300 kms ') and the ICL (with VLOs.b - 1000 
kms ' and cr\, ^ 900 kms '). Furthermore, in Sect. |5] we found 
that the halo and IC components were characterised by specihc 


PN numbers (a parameters) that differed by a factor of three, and 
by different shapes of their PNLFs. 

These results demonstrate the coexistence of two distinct PN 
progenitor stellar populations in this region of the Virgo clus¬ 
ter core: the M87 halo and the ICL. These two populations have 
very different surface density distributions, the M87 halo is de¬ 
scribed by an n = 11.8 Sersic law, while the ICL follows a shal¬ 
low power law oc Rr with y in the range [-0.34,-0.04]. We 
also have external information on the metallicities and ages of 
both components. At B ~ 35kpc, the mean metallicity of the 
M87 halo obtained with population synthesis models from multi¬ 
colour photometry is ~ 0.7 solar, with a shallow outward gradi¬ 
ent, a nd the mean age is ~ lOGyr (iLiu et al.ll2005t [Montes et ^ 
l2014[l . On the other hand, the metallicity and age distributions of 
ICL red giants in a held atB ~ 190 kpc from HST ACS star pho- 
tomet ry are dominated by m etal-poor ([M/H1< -1), > lOGyrold 
stars ([Williams et al.lKOO^ . Because of the large velocities and 
shallow surface density prohle of the IC stars, these IC popula¬ 
tion parameters are likely to be similar in the radial range probed 
by our observations, B ~ 50 - 140 kpc, whereas the M87 halo 
stars might reach ~ 0.5 solar in the outer regions if the outward 


_ regii _ 

gradient continues, as inferred bv iLiu et’aD (l2005h . 

Currently there is no good theoretical understanding of how 
the properties of a PN population are related to the metal¬ 
licity and age of a stellar population. Observationally, star¬ 
forming and bulge populations have a numbers such as those 
we hnd for the ICL, while only the most massive early-type 
galaxi es have a numbers as low as we hnd for the M87 
halo (iBuzzoni et al.l [20061 : ICortesi et alJ l2013l) . The primary 
driver is believed to be increased mass loss at high metallic¬ 
ities. The PNLFs are empirically found to steepen from star¬ 
forming to old metal-r ich populations ([Ciardullo et al.l 120041 : 
iLongotordi et al.l 1201 3ll . The PN LFs of Local Group star ¬ 
formin g galaxies, such as th e SMC ( Jacobv & De Marco[|2002l) . 
LMC (IReid& Par kdiZOlOi). M33 (ICiardullo et al.l l2004 . md 
NGC 6822 (iHernandez-Martmez & Pefia 20091) . furthermore 
show a ‘dip’ 2-2.5 mag down from the PNLF cut-off for the 
LMC, M33, and NGC 6822, and 4 mag down from the cut¬ 
off for the SMC. A tentative model for this feature is the su¬ 
perposition of a faint PN population with a brighter popula- 
tion of more massive cores from a younger stellar population 
(iRodrfguez-Gonzalez et alJl2014l) . We can speculate that as the 
brighter population fades in older and/or more metal-rich pop¬ 
ulations, the dip might move towards brighter magnitudes. This 
could explain why in the Virgo IC population we find the dip 1 - 
1.5 mag down from the cutoff. No other PN population with this 
PNLF is known; however, PNLFs as deep as for M87 have only 
been obtained in the Local Group so far. Clearly, more observa¬ 
tional and theoretical work on the nature and location of the dip 
in the PNLF is needed. 


6.2. The ICL in Virgo: nature and number of its progenitor 
gaiaxies 

The combined properties of the Virgo ICPN population, which 
include the fairly small inferred bolometric luminosity, the rel¬ 
atively large a parameter, and the dip in the PNLF as well as 
the low mean metallicity from [Williams et ^ ([20071) . appear to 
be most readily explained if this population derives from a faded 
population of low-luminosity, low-metallicity, star-forming or ir¬ 
regular galaxies, such as M33 or the LMC, which are very dif¬ 
ferent from M87 itself. 

In Section 15.11 we determined the total V-band and bolo¬ 
metric luminosities of the IC component sampled by our survey 
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fields: Licl = 0.53 x 10^°Lo, and LicL.boi = 1.1 x 10'°Lo,boi- 
Using the total V-band luminosities for M33 and the LMC listed 
in NED, Lm33=3.65x10^Lo and Llmc=1-26x10^Lo 5 we find that 


the IC stars sampled in our survey fields correspond to ~ 1.5 
M33-like galaxies or ~ 4 LMC-like galaxiefl 

We can now also check whether the Virgo ICL associated 
with the ICPN population could be related to the blue GC pop¬ 
ulation that is found around M87, that has a shallower and 


trace the stellar halo light (Cote et al. 

2001; Tamura et alj|2006l: 

IStrader et al.ll201ll iForte et alJ2012: 

Durrell et al. 20141). To do 


this, we need to estimate the total number of GCs as sociated 
with M33 and LMC-like systems. iHarris et aP (1201 3h studied 
GC populations in a large sample of galaxies and analysed the 
correlation of the total number of GCs, Nqc, with global galaxy 
properties and type. They find that Ngc increases roughly in 
direct proportion to host galaxy luminosity, with a scatter of 
a factor of ~ 2.5. For an LMC-like system with luminosity 
Llmc=1-3x10^Lo the expected mean number of GCs is Noe ~ 
20, while for an M33-like galaxy with Lm33=3.7x10®Lo the ex¬ 
pected mean number of GCs is Nqc ~ 60. This leads to an esti¬ 
mated number of blue GCs associated with the sampled ICL of 
IVgc.icl ~ 80-90, with a scatter of a factor ~ 2.5. If a fraction of 
the ICL is due to the accretion of even lower lumi nosity galax¬ 
ies, the estimated number of GCs would increase (IHarris et alJ 
l2013HCoccato etalJl2013h . 


The recent s urvey of the GC pop ulation in the Virgo cluster 
around M87 by iDurrell et al.l (1201 4l) showed the presence of an 
ICGC population, mostly associated with blue GCs (see further 
discussion in Sect. 16.31) . This intracluster component has a den¬ 
sity equal to l,bGC,tot - 0.2^g ggarcmin“^. In our surveyed region, 
this would lead to a total number of 100-430 ICGCs. This is 
larger than but consistent within the uncertainties with the value 
estimated above, suggesting that a substantial fraction of the blue 
GC population around M87 could have been accreted with the 
galaxies that we now see in the ICL. 


6.3. Is there an iniraclusler component of globular clusters 
around M87? 

There has been some controversy in the recent literature about 
the existence of intracluster GCs in the halo of M87. The most 
extensive photometric study of the distri bution of GCs in the 
Virgo cluster so far was carried out by iDurrell et al.l (1201 4l) 
as part of the NGVS. They studied density maps of the GC 
population, selected using colour criteria, and statistically ac¬ 
counted for the contamination to the GC sample by subtract¬ 
ing a modelled map for the expected ba ckp ound, from both 
Milky Way stars and background galaxieo The blue GCs in 
their map h ave a shallower and more extended profile than 
the red GCs. IDurrell et al.l (1201 4h also found that the total GC 
(blue plus red) density profile is in good agreement with the 
number density profile of photometrically selected PNs from 
iLongobardi et al.l (1201 3l) . including a change of slope and a flat¬ 
ter profile at large radii. They suggested that their blue GCs at 
distances > 215 kpc are part of t he intracluster component of 
Virgo. Cosmological simulations (iDolag et al.l 1201 Ol; ICui et aTI 


* The total luminosity at all radii corresponding to the detected IC 
stars is much larger; their large measured velocity dispersion implies 
that the orbits of these IC stars reach to much larger radii in the cluster. 

We note here that the M87 PN sample cannot be contaminated by 
Milky Way halo PNs, as these would have [OIII]/i5007A fluxes about 
12 mag brighter than M87 PNs. 


1201 4t ICooner et al.ll2014ll predict that the density of this com- 
po nent would then incre ase inwards. 

IDurrell et al.l (1201 4l) also found evidence for a spatial asym¬ 
metry of GCs surrounding M87 for major axis distances larger 
than 20', with an excess of tracers in the NW region (mostly the 
blue population). In Sect. |4]we studied the distribution of M87 
halo and IC PNs separately for the NW and SE. We find no clear 
evidence of asymmetry in either the halo and ICL within major 
axis distance ~ 20' (Fig. [loll. Inside this radius, both PN and GC 
number density profiles are consistent with a symmetric halo and 
IC distribution. For the halo component, this result is significant, 
given the number of tracers and radial extent, and indicates, that 
if the halo was subject to accretion events, these were not re¬ 
cent. For the ICL, we may expect asymmetries, given the longer 
timescales involved in IC accretion events, but we may not have 
a l arge enough sample o f ICPNs to see them. 

IStrader et akl (1201 ll) carried out a spectroscopic study of the 
GCs around M87, using colour criteria to select their candidates. 
In the same colour and magnitude range populated by globular 
clusters, 0.55 < (g' - /') < 1.15 and 20 < g' < 24, there is 
however a large contribution from foregro u nd Mi lky Way halo 
stars. To mitigate this effect, IStrader et al.l (1201 ll) used a com¬ 
bination of photometry, radial velocity, and HST imaging in¬ 
formation. However, they considered all objects with velocity 
Vlos < 150 kms"' to be stars. In the most ambiguous range, 
i.e. 150kms“’ < Vlos < 350kms“', they classified all uncertain 
objects, for which a clean separation between contaminants and 
GCs could not have been done, as stars. Based on the remain¬ 
ing sample, they reported that the number density profile of the 
spectroscopically confirmed GCs showed no evidence for a tran¬ 
sition between a halo and IC component, either as a sharp trun¬ 
cation of the halo, or a flattening of the GC number density pro¬ 
file at large radii. From the sample kinematics, they observed 
that their GCs around M87 have velocity dispersion in the range 
300 < cr < 500 kms * out to 190 kpc, with ~ 500 kms ' for 
the GCs population at 190 kpc significantl y smaller than the ve¬ 
locity dispersion of V irgo cluster galaxies (iBinggeli et al.ll 19931: 
IConselice et al.ll20()lh . However, from the PN phase-space dis¬ 
tribution (VLos,PN>f^PN) in Fig- in we see that a large fraction of 
IC stars near M87 have velocities Vlos < 350 kms '. This sug- 
gests that the lack o f evidence for the IC component reported by 
IStrader et al.l (1201 ll) could be caused by the velocity threshold 
Vmin = 350 kms“* imposed on the GC sample, which is needed 
to prevent the contamination from Milky Way halo stars, but may 
also remove many of the ICGCs from their analysis. 

6.4. Relation between BCG and ICL 

When studying central galaxies in clusters, one of the main ques¬ 
tion is to establish where the ICL begins and where the associ¬ 
ated BCG ends, or whether any distinction is to be made at all. 

For M87, the differences in the density profiles and veloc¬ 
ity distributions of the halo and ICPN populations, as well as in 
their a-parameters and PNLFs, are sufficient to argue that the 
two components must be considered to be separate stellar popu¬ 
lations, with different metallicities and star formation histories, 
and not as a continuum. As discussed above, published stellar 
population data suggest that the halo stars are older and more 
metal-rich than the ICL (see Sects. |5]|6T]for more details). 

In more distant BCGs where a kinematic decomposition be¬ 
tween BCG halo and ICL is not available, the presence of an ad¬ 
ditional dynamical component in BCGs is usually i nferred from 
a change of slope at large radii in the SB profile (IZibetti et al.l 
I 2 OO 5 I : iGonzalez et al.lI 2 OTO ID’ Souza et alJ 120141) . Photometric 
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properties or colours are obtained by treating the two compo¬ 
nents as a continuum because no differentiation between the un¬ 
derlying stellar populations is normally possible. 

Using a particle tag ging method t o ana lyse galaxy clusters 
in ACDM simulations, ICooner et al.l (1201 4h consider the BCG 
and ICL as a single entity consisting of all stars, which are not 
bound to any cluster subhalos. They then split the BCG stars 
into accreted stars and in situ stars, and find that the large ma¬ 
jority of BCG stars are accreted stars. They find double-Sersic 
surface density profiles in their simulated BCGs, where the in¬ 
ner component {R < 200kpc) is dominated by ‘relaxed’ accreted 
components, and the outer comp onent by ‘unrelaxed’ accreted 
components. ICooner et al.l (|2014l) argue that the accreted/in situ 
separation is physically meaningful and that the ICL should nat- 
uraly be considered as a continuation of the BCG to low sur¬ 
face brightness because both components are formed by similar 
mechanisms. 

In contrast, a dynamical approach based on the velocity dis¬ 
tributions of diff use light particles i n hydrodynami c cosmolog¬ 
ical simulations dPolag et alj|2010l: ICui et al.fl2Q14ll is found to 
separate these stars into two components: one bound to the clus¬ 
ter potential and the other bound to the BCG. The resulting BCG 
and diffuse ICL are formed on different timescales, and the sim¬ 
ulated stars associated with the two components are d ifferent in 
terms of spatial distribution, ages, and metallicities. ICui et alJ 
(l2Q14ll showed that it is possible to dynamically differentiate 
between halo and IC particles using the particles’ binding en¬ 
ergies. Stars with high binding energy, which end up belong¬ 
ing to the BCG, were subjected to relaxation and merging pro¬ 
cesses such that the gravitational potential changed so quickly 
that these stars lost memory of the kinemat ics of their progeni¬ 
tors dMurante et al.ll200^lDolag et al.ll2M^ . On the other hand, 
stars with lower binding energy, that belong to the diffuse com- 
ponent, still reflect the dynamics of the satellite galaxies. Both 
iDolag et^ d201f)l) and ICui et al.) d2014l) also observed that the 
slope of the surface brightness profile associated with the two 
components change, with the halo profile being steeper than the 
ICL profile. 

It is likely that the distinct BCG and IC components found 
in the hydrodynamical simulations are related to the relaxed and 
unrelaxed accreted components in the particle tagging analysis, 
but the inclusion of baryonic processes in the former may ac¬ 
centuate the differences found between BCG and I CL. If we 
associ ate the BCG and ICL of iDolag et^ d2010ll : ICui et al.l 
( 20141) with the rela xed and unrelaxed accreted components of 
ICooner et al.l d2014l) . the progenitors of the stars in the steeper 
Sersic (relaxed) component would be accreted from more mas¬ 
sive systems at higher redshifts, while the stars in the shal¬ 
lower and more extended ICL (unrelaxed) component would 
come from the accretion of less massive systems at lower red- 
shifts. More massive progenitors domin ate the diffuse light 
in simulated clusters close to the centre dMurante et alJl2()TO 
iPuchwein et al.ll20T^ because they move further inwards by dy¬ 
namical friction. They cause stronger relaxation of the gravita¬ 
tional potential, and if accreted early, they have more time to 
relax. 

To summarise, recent simulations show that a distinction can 
be made between stars that trace the cluster potential and stars 
bound to the BCG, based on the physical properties and binding 
energies of the accreted progenitors. From the study of the PN 
population around M87, we have shown the coexistence of two 
discrete components in the Virgo cluster core, tracing different 
stellar populations, in agreement with these predictions. While 
the PN population for the IC component around M87 indicates 


low-mass dwarf and star-forming galaxy progenitors, the stellar 
halo has higher metallicity, ~ 0.7 solar, indicating more massive 
progenitors. This bimodality in the progenitors may be the root 
of the bimodality in the kinematics and density profiles of the 
M87 halo and the ICL. 

However, we note that this kind of bimodality need not occur 
in every cluster of galaxies. For example, it is plausible that, for 
a more continuous distribution of progenitor masses and a more 
uniform distribution of binding energies of the debris stars, the 
final BCG plus ICL system would show continuous radial gradi¬ 
ents in kinematics and stellar population properties, rather than 
appear as the sum of several discrete components. It is possible 
that NGC 6166 in the Abell 2199 cluster is closer to this situa¬ 
tion: the velocity dispersion in the high surface brightness halo 
of this BCG was recently measured to increase up to the clus¬ 
ter vel ocity dispersion of ~ 800km.s ^ at 100” from the galaxy 
centre (iBender et alJ (|2014|) : see also lKelson et al.l (l2002h l. 


7. Summary and conclusions 

We obtained spectra for 287 PNs in the outer regions of the 
nearby elliptical galaxy M87, of which 211 are located be¬ 
tween distances 40 kpc to 150 kpc from the galaxy centre. 
Spectra were acquired with the FLAMES spectrograph in the 
GIRAFFE-i-MEDUSA configuration, with spectral resolution of 
R = 22500. We observed 14 different ELAMES plate con- 
figurations, using candid ates from the catalogue described in 
iLongobardi et al.l (l2013l) . The spectroscopic survey aimed at 
measuring the LOS velocities of PNs in the transition region be¬ 
tween the galaxy’s stellar halo and the ICL. We identified PNs 
through their narrow and symmetric, redshifted [OIII]/15007A 
emission line, with no or negligible continuum, and verified with 
the second [OIII]/14959A emission line. Spectra were measured 
for PNs in the magnitude range from m5oo7=26.3 down to 28.8. 
This is the largest spectroscopic sample of PNs at such galactic 
radii for a central galaxy, in the number of tracers and magnitude 
depth. 

The area covered by the survey allowed us to trace the tran¬ 
sition between the M87 halo and ICL in the Virgo cluster core. 
The coexistence of these two components is shown by the bi¬ 
modality of the LOSVD, whose strong asymmetric wings make 
it deviate from the near-Gaussian LOSVD typical of early-type 
galaxies. We separated halo and ICPNs by studying the projected 
phase-space distribution. We implemented a robust technique to 
measure the velocity dispersion of the M87 halo, separating its 
velocity distribution from the broader component, the ICL. We 
identified 243 PNs for the M87 halo and 44 ICPNs. We found 
that the logarithmic number dens ity profile for the halo P Ns fol¬ 
lows the V-band SB profile from iKormendv et al.l (l2009l) . while 
the IC number density profile decreases towards large radii as a 
power-law /icl R^ with y in the range [-0.34, -0.04]. 

The total PN surface density profile at large radii is flatter 
than the surface brightnes s profile because of the presence of 
the IC component (see also ILongobardi et al.ll2013h . which con¬ 
tributes ~3 times more PNs per unit luminosity than the halo 
population. We find luminosity-specific PN numbers ahaio = 
(1.06 + 0.12) X 10-^ NpNLgboi Q-icL = (2.72 + 0.63) x lO^* 
NpnLq*{,q[ for the M87 halo and IC PN population, respectively. 
This is consistent with the known existence of a gradient towards 
bluer colours at large radii because of the increased contribution 
of ICL at large distances and its lower metallicity compared to 
the halo population. 
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The spectroscopically confirmed PNLFs for both the halo 
and IC PNs have a steeper slope towards fa int magnitudes 
than i s predicted by the analytical formula of ICiardullo et alJ 
( 19891). confirming the result from the photometric sample 
(iLongobardi et al.ll2013h . This steepening is consistent with an 
old stellar population dominated by PNs with low-mass cores. 
The PNLF of the ICPN population has a slightly shallower gradi¬ 
ent than the M87 halo PNLF, and in addition shows dip at about 
~ 1 - 1.5 magnitudes from the bright cut-off. This dip is an evo¬ 
lutionary feature observed in star-forming systems, such as M33 
and the Magellanic clouds, and may be related to rapidly evolv¬ 
ing PNs with massive central cores. The presence of the dip in 
the ICPNLF, but not in the M87 halo PNLF, provides additional 
evidence for intrinsic differences between the halo and IC parent 
stellar populations. 

Using PNs as tracers we showed that the stellar halo of the 
BCG galaxy M87 is distinct from the surrounding ICL in its 
kinematics, density profile, and parent stellar population, consis¬ 
tent with the halo of M87 being redder and more metal-rich than 
the ICL. We note that the ICL in our surveyed fields corresponds 
to about four times the luminosity of the LMC, spread out over 
a region of ~ lOOkpc diameter. It is remarkable that population 
properties can be observed for such a diffuse component. 

In the Virgo cluster, BCG halo and ICL cannot be consid¬ 
ered as components with a gradual transition in their kinematics. 
This supports results from analysis of galaxy cluster simulations, 
which suggest that the IC component in Virgo consists of un¬ 
relaxed accreted stars bound to the cluster potential, while the 
stellar halo of M87 appears to be described as a relaxed accreted 
component bound to the galaxy itself. Based on its PN popu¬ 
lation properties, we propose that the progenitors of the Virgo 
ICL were low-mass, star-forming galaxies, which may also have 
brought with them a significant fraction of the blue GC popula¬ 
tion seen in the outer regions of M87. 
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